{ "cells": [ { "cell_type": "markdown", "id": "586b41e3", "metadata": {}, "source": [ "# 04 - Raven Discharge Ensemble Forecasts" ] }, { "cell_type": "markdown", "id": "85360e9f", "metadata": {}, "source": [ "## Learning goals of this module\n", "- Learn how to do a basic analysis, comparing observed discharge with forecast discharge, produced by Raven and forced by GEPS, GEFS and IFS.\n", "- Learn about ensemble forecast verification.\n", "- Learn how to calculate and visualize simple error metrics, such as the __rank_histogram__ and the __continuous_ranked_probability_score__.\n", "\n", "## Assumptions\n", "- We assume you are familiar with the concept of a __probabilistic__ forecast. \n", "- We assume you are familiar with the concept of a __hydrological__ model.\n", "\n", "### Reference to ensemble forecast products\n", "\n", "- [Global Ensemble Prediction System (GEPS)](https://open.canada.ca/data/en/dataset/6d9dd2f8-202e-58cb-a110-e2168832aacb)\n", "\n", "- [Global Ensemble Forecast System (GEFS)](https://www.ncei.noaa.gov/products/weather-climate-models/global-ensemble-forecast)\n", "\n", "- [Integrated Forecasting System (IFS)](https://www.ecmwf.int/en/forecasts/documentation-and-support/changes-ecmwf-model)\n", "\n", "### Reference to Raven\n", "[Raven Hydrological Model](https://raven.uwaterloo.ca/)\n" ] }, { "cell_type": "markdown", "id": "7a6aac68", "metadata": {}, "source": [ "## Run imports and set-up logging" ] }, { "cell_type": "code", "execution_count": 1, "id": "edc0b447", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2026-06-18 22:49:17,077 - INFO - Running Veriflow version 0.1.0\n" ] } ], "source": [ "import logging\n", "import sys\n", "import warnings\n", "from pathlib import Path\n", "\n", "from dotenv import load_dotenv\n", "\n", "from veriflow import run_pipeline\n", "from veriflow.constants import VERSION\n", "\n", "# add project root (parent of notebook folder) to path\n", "sys.path.append(str(Path(\"..\").resolve()))\n", "\n", "\n", "from verification_plots import (\n", " crps_plot,\n", " forecast_timeseries_plot,\n", " rank_histogram_3d_plot,\n", " rank_histogram_plot,\n", " scatter_plot,\n", ")\n", "\n", "# Reload automatically\n", "%load_ext autoreload\n", "%autoreload 2\n", "\n", "\n", "warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n", "warnings.filterwarnings(\"ignore\", category=FutureWarning)\n", "\n", "load_dotenv(dotenv_path=\"tutorial.env\", override=True)\n", "\n", "base_config = Path(\"config\")\n", "base_config.exists()\n", "\n", "logging.basicConfig(\n", " level=logging.INFO,\n", " format=\"%(asctime)s - %(levelname)s - %(message)s\",\n", " handlers=[logging.StreamHandler()],\n", ")\n", "logging.info(f\"Running Veriflow version {VERSION}\")" ] }, { "cell_type": "markdown", "id": "5dd248c2", "metadata": {}, "source": [ "## Inspecting the _veriflow_ pipeline configuration\n", "1. Open the config file in the \"config\" directory. The name of the file is identical to the name of the notebook.\n", "2. Inspect each of the sections to gain an understanding of what this configuration is about." ] }, { "cell_type": "markdown", "id": "ec89e13b", "metadata": {}, "source": [ "## Running the _veriflow_ pipeline" ] }, { "cell_type": "code", "execution_count": 2, "id": "b339ac88", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2026-06-18 22:49:17,331 - INFO - Successfully initialized the configuration. \n", "\t verification_period_start = 2025-05-20 00:00:00 \n", "\t verification_period_end = 2026-06-01 00:00:00\n", "2026-06-18 22:49:17,340 - INFO - Start getting data from FewsWebservice.\n", "2026-06-18 22:49:18,007 - INFO - Successfully got observed data from FewsWebservice.\n", "2026-06-18 22:49:18,008 - INFO - Start getting data from FewsWebservice.\n", "2026-06-18 22:49:18,693 - INFO - Successfully got simulated_GEPS data from FewsWebservice.\n", "2026-06-18 22:49:18,694 - INFO - Start getting data from FewsWebservice.\n", "2026-06-18 22:49:19,598 - INFO - Successfully got simulated_GEFS data from FewsWebservice.\n", "2026-06-18 22:49:19,600 - INFO - Start getting data from FewsWebservice.\n", "2026-06-18 22:49:20,632 - INFO - Successfully got simulated_IFS data from FewsWebservice.\n", "2026-06-18 22:49:20,637 - INFO - Successfully loaded all data from sources.\n", "2026-06-18 22:49:20,703 - INFO - Successfully computed CrpsForEnsemble for verification pair QR_GEPS.\n", "2026-06-18 22:49:20,774 - INFO - Successfully computed CrpsForEnsemble for verification pair QR_GEFS.\n", "2026-06-18 22:49:20,852 - INFO - Successfully computed CrpsForEnsemble for verification pair QR_IFS.\n", "2026-06-18 22:49:20,949 - INFO - Successfully computed RankHistogram for verification pair QR_GEPS.\n", "2026-06-18 22:49:20,991 - INFO - Successfully computed RankHistogram for verification pair QR_GEFS.\n", "2026-06-18 22:49:21,038 - INFO - Successfully computed RankHistogram for verification pair QR_IFS.\n", "2026-06-18 22:49:21,040 - INFO - Verification pipeline completed successfully.\n" ] } ], "source": [ "ods = run_pipeline((base_config / \"04_raven_elbow_discharge_ensemble.yaml\", \"yaml\"))" ] }, { "cell_type": "code", "execution_count": 3, "id": "4c458832", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "OutputDataset with 3 verification pairs: [VerificationPair(id='QR_GEPS', obs='observed', sim='simulated_GEPS', variable='discharge'), VerificationPair(id='QR_GEFS', obs='observed', sim='simulated_GEFS', variable='discharge'), VerificationPair(id='QR_IFS', obs='observed', sim='simulated_IFS', variable='discharge')]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ods" ] }, { "cell_type": "markdown", "id": "0d7caa3a", "metadata": {}, "source": [ "## Evaluating the results in the _veriflow_ `OutputDataset`" ] }, { "cell_type": "markdown", "id": "cb8dd8d5", "metadata": {}, "source": [ "Verification metrics and results can contain a level of abstraction. Although these abstractions can reveal important information about forecast quality, a basic \"eyeball verification\" is often the best and intuitive way to start your verification exercise. You'll likely find strengths and weaknesses in your forecasts early on, without directly diving into levels of abstraction. In addition, a solid visual inspection may help you later on in understanding or explaining the more abstract results." ] }, { "cell_type": "markdown", "id": "d178a47e", "metadata": {}, "source": [ "### 1 - Visual inspection of observed and forecast data\n", "A good starting point for \"eyeball\" verification is simple: just looking at your observations and forecasts in a visual way. Use the interactive elements in the plots below to zoom, pan and compare the results of our 3 NWP products." ] }, { "cell_type": "code", "execution_count": 4, "id": "ab797cc8", "metadata": {}, "outputs": [], "source": [ "stations = ods.get(ods.verification_pairs[0]).coords[\"station\"].values\n", "lead_times = ods.get(ods.verification_pairs[0]).coords[\"lead_time\"].values" ] }, { "cell_type": "code", "execution_count": 5, "id": "9b157dad", "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "forecast_timeseries_plot(ods)" ] }, { "cell_type": "markdown", "id": "546ab3d5", "metadata": {}, "source": [ "### 2 - Visual inspection with a scatter plot per lead time\n", "Another great tool for \"eyeball\" verification is the scatter plot. The scatter plot is relatively easy to understand, but is slightly more abstract than the visualization above. \n", "\n", "For all forecasts in our `OutputDataset`, we collect all realizations of the ensemble at a specific lead time. You can use the `lead_times` variable (a list of `np.timedelta64` instances) to slice the data at a specific lead time. For that given slice, we can make the scatter plot.\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "33bbf33f", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "scatter_plot(ods, lead_time=lead_times[1])" ] }, { "cell_type": "markdown", "id": "e76c0fb8", "metadata": {}, "source": [ "### 3 - Looking into the Continuous Ranked Probability Score\n", "\n", "Next, we'll look into the results of the continuous ranked probability score. Before you continue to the next section, we'll look into the documentation and do a tutorial on the CRPS for ensemble forecasts. The _veriflow_ package relies on _scores_ (developed by the Bureau of Meteorology, Australia) for computation of various scores. \n", "\n", "1. Read the [API documentation](https://scores.readthedocs.io/en/stable/api.html#scores.probability.crps_for_ensemble) the __crps_for_ensemble__ function, which is used under the hood in _veriflow_. \n", "2. Run through the [scores tutorial](https://scores.readthedocs.io/en/stable/tutorials/CRPS_for_Ensembles.html) on the __crps_for_ensemble__ function. You can run it in Binder (link on top of the page), or view the static view.\n", "\n", "\n", "- Q1: what attribute(s) of forecast quality can be measured by the CRPS?\n", "\n", "
\n", "Show suggested answers for Q1\n", "The CRPS captures multiple attributes: accuracy, reliability and sharpness. Can you reason why?\n", "\n", "
\n", "\n", "- Q2: what is the best possible CRPS score?\n", "\n", "
\n", "Show suggested answers for Q2\n", "The best possible outcome of CRPS is 0. In this case the forecast has maximum sharpness: all ensemble members exactly predict the observed outcome.\n", "
\n", "\n", "\n", "- Q3: computing the CRPS over just one realization (i.e. a deterministic forecast) yields the exact same result as computing the ... for a deterministic score?\n", "\n", "
\n", "Show suggested answers for Q3\n", "The absolute error. The CRPS is a probabilistic generalization of the absolute error. When taking the mean of the CRPS over all forecasts, it is equal to the mean absolute error when the number of realizations is 1.\n", "
" ] }, { "cell_type": "code", "execution_count": 7, "id": "41bda300", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "crps_plot(ods)" ] }, { "cell_type": "markdown", "id": "c777af8a", "metadata": {}, "source": [ "### 4 - Looking into the Rank Histogram\n", "\n", "Next, we'll look into the results of the rank histogram. Before you continue to the next section, we'll look into the documentation and do a tutorial on the rank histogram for ensemble forecasts. The _veriflow_ package relies on _scores_ (developed by the Bureau of Meteorology, Australia) for computation of various scores. \n", "\n", "1. Read the [API documentation](https://scores.readthedocs.io/en/stable/api.html#scores.probability.rank_histogram) the __rank_histogram__ function, which is used under the hood in _veriflow_. \n", "2. Run through the [scores tutorial](https://scores.readthedocs.io/en/stable/tutorials/Rank_Histogram.html) on the __rank_histogram__ function. You can run it in Binder (link on top of the page), or view the static view.\n", "\n", "- Q1: what attribute(s) of forecast quality can be measured by the Rank Histogram?\n", "\n", "
\n", "Show suggested answers for Q1\n", "The rank histogram primarily measures reliability, although bias in forecasts will show up in the rank histogram as well.\n", "
\n", "\n", "- Q2: what does a perfect Rank Histogram look like?\n", "\n", "
\n", "Show suggested answers for Q2\n", "The rank histogram shows a uniform (flat) distribution, indicating forecast and observed frequencies are equal.\n", "
\n", "\n", "- Q3: does a perfectly uniform/flat rank histogram always correspond to a good forecast? Can you come up with a hypothetical forecast that has a perfect rank histogram, but is still bad?\n", "\n", "
\n", "Show suggested answers for Q3\n", "No, although a perfect rank histogram indicates statistical reliability, the forecast may still have poor resolution (the ability to discriminate between situations that have different observed outcomes).\n", "\n", "Imagine an ensemble forecast that ignores current conditions and simply samples from the historical distribution of streamflow for that day of year.\n", "\n", "For example, every day in June the ensemble consists of random draws from the last 30 years of June observations.\n", "\n", "This forecast can also produce a nearly uniform rank histogram because the observations come from the same climatological distribution. However, it has zero ability to distinguish today's conditions from any other June day.\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "72b2b8af", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rank_histogram_plot(ods, station=stations[0], lead_time=lead_times[4])" ] }, { "cell_type": "markdown", "id": "d283f52c", "metadata": {}, "source": [ "### 5 - Looking into the Rank Histogram in 3D (advanced)\n", "This might seem relatively abstract or advanced, but the 3D rank histogram shows the rank histogram for all lead times. In the current notebook, the number of samples is too low to make a proper analysis of the results, but with many samples this can be an interesting and useful visualization.\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "fb7d77ba", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rank_histogram_3d_plot(ods, station=stations[0])" ] } ], "metadata": { "kernelspec": { "display_name": "veriflow", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.4" } }, "nbformat": 4, "nbformat_minor": 5 }